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We consider recently proposed higher order gravity models where the action is built from the 
£S) , Einstein-Hilbert action plus a function f(G) of the Gauss-Bonnet invariant. The models were 

previously shown to pass physical acceptability conditions as well as solar system tests. In this 
O |' paper, we compare the models to combined data sets of supernovae, baryon acoustic oscillations, 

\ and constraints from the CMB surface of last scattering. We find that the models provide fits to the 

data that are close to those of the LCDM concordance model. The results provide a pool of higher 
order gravity models that pass these tests and need to be compared to constraints from large scale 
structure and full CMB analysis. 
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Higher order gravity models have been proposed among possible causes of late-time cosmic acceleration [l[ . These 
Oh models are built out of higher order curvature invariants that yield generalized field equations with a coupling between 
the mass content of the universe and the space-time curvature that produces a late-time self accelerating phase. A 
large body of papers have been devoted to the so-called f(R) models while a smaller fraction study models based 
' on invariants built out of the Ricci and Riemann tensors Q. 

In addition to the phenomenology of an accelerating cosmic expansion, higher-order gravity models have theoretical 
motivations within unification theories of fundamental interactions, and through field quantization on curved space- 
, times 0j H- Indeed, higher-order invariants appear automatically in most quantum gravity proposals 043]' string 
J> ■ theories [all], supergravity [H, |Tl|, and loop quantum gravity/cosmology [l2rll3 ]. In these theories, high-order loop 
ON ! corrections on curved space-time are related to higher-order combinations of the Riemann curvature invariants. Higher- 
order invariants have been actively discussed in relation to avoiding cosmological curvature singularities (see 15[ and 
references therein). In some of these theories, the invariants are regrouped in a topological invariant combination 
called the Gauss-Bonnet invariant, denoted as G. This combination gives a theory free of unphysical states [H, EH El]] • 
In this paper, we study cosmological constraints on some models where the action is made of the Einstein-Hilbert 
action plus a function f(G) of the Gauss-Bonnet invariant. We focus on models that have been shown in previ- 



ous literature to be free of ghost instabilities and superluminal propagations [181 ] in cosmological homogeneous and 
isotropic backgrounds. We perform coordinate transformations in order to write the dynamical equations in a form 
integrable using numerical schemes and then compare the models to recent observations of supernova magnitude- 
redshift data, distance to the CMB surface, and Baryon Acoustic Oscillations (BAO), including Hubble Key project 
and age constraints. 

2. /(G) MODELS 

The models that we investigate in this paper are derived from varying the action 
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with respect to the metric 5g a p, where 

G = R 2 - 4R al3 R al3 + R aMS R aM5 (2) 

is the Gauss-Bonnet invariant, R is the Ricci scalar, R a p is the Ricci tensor, R a /3-ys is the Riemann tensor, L m and 
Lrad are the matter and radiation energy Lagrangians, respectively. We will work in units with reduced Planck mass 
M 2 { = (SttGn)^ 1 = 1. The corresponding field equations read 

ift)]V 7 V*/ G 

+ (G/g- f)g a p + R a p -^g a /3R = T al3 , (3) 

where we use the definition fa = 

Now using the Friedmann-Lemaitre- Robertson- Walker flat metric 

ds 2 = -di 2 + a{tfdx 2 (4) 

and Universe filled with matter and radiation, one derives the generalized Friedmann equation 

3H 2 = Gf G - fG - 2AH 3 f G + Pm + Prad . (5) 

where p m and p rac [ are the matter and radiation energy densities, respectively, H = a/a and a dot represents d/dt. 
We also note that in terms of H, 

R = 6{H + 2H 2 ) (6) 

and 

G = 24H 2 (H + H 2 ). (7) 

We explore explicit models in the next sections. 

3. RECENTLY PROPOSED VIABLE /(G) MODELS 

3.1. f(G) models proposed by De Felice and Tsujikawa 

In Ref. [20| . the authors imposed certain conditions on the function /(G) and its derivatives. The most important 
condition being d 2 f/dG 2 > in order to ensure the stability of a late-time de Sitter solution as well as the existence 
of a radiation/matter dominated epochs preceding a late-time accelerating phase. Other additional regularity and 
viability conditions in [20( single out the following /(G) functions 



Model DFT-A: /(G) = A^= arctan (—) - -^/Chln (l + -J)- a\^fG~ , (8) 

yGo ^Go/ 2 \ Gq/ 

Model DFT-B: /(G) = \-== arctan (—) - a\^/C^ : (9) 

VGo ^Go^ 

Model DFT-C: /(G) = \VG\n [cosh (— )] - a\^/G~ a , (10) 

L VGo 

where A, Go and n are positive constants and a is a constant. These functions were derived from the integration of the 
following second order derivatives that satisfy the condition J.gg > for all values of G along with other regularity 
conditions [20(. 

Model DFT-A: J.gg = -™ — , n = 1 (11) 

Gl' 2 [l + {G*/GIY\ 

2A 

Model DFT-B: /, GG = — , n = 2 (12) 

Gl' 2 [l + (GVGg)»] 

Model DFT-C: f,GG = -375 [1 - tanh 2 (G/G )], (13) 
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Varying the corresponding actions with respect to the metric, the generalized Friedmann equations follow as: 
Model DFT-A: 



3H 2 



G (1152H 6 (H{4:H' + 3H + H") - a(H + H 



/\2\ 
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i n 
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(14) 



Model DFT-B: 



3H 2 = {aGl - 576G 2 H 6 (H / + H)(5H' + H) + 1152aG 2 H 6 (H' + H) 2 

+331776tf 12 (ff' + H)\l + a)) ^G~ X / (G 2 + G 2 ) 2 + + 3H &™ d + A 



(15) 



Model DFT-C: 



G 



3H 2 = A 24 sinh — G H 3 cosh ( — ){H' + H) 



Go 

-576# 6 (3£T' 2 - H(4H' + H")) 
G 



G 



Go 



Go cosh ( — I ( a - In ( cosh ( — 



Go 



G 



Go 



3-//q £l<m 



3Hq £l ra d 
Tan 



A. 



(16) 



where A = a\\/G~o and we define N — In a with ' = d/dN. As discussed in [19|, |21], these higher order gravity models 
have generalized Friedmann equations that are very stiff ordinary differential equations (ODEs) that require us to 
use well-adapted numerical codes and logarithmic variables. Thus, we replaced H in the ODEs by the logarithmic 
variable u = \n(H/(l) where /t is a constant normalizing the Hubble parameter. This allows one to write the source 
terms, see Bill, as 



u = ln (Q rad e- 4N +u m e- 3N )/2 



(17) 



and 



8irG po 



m ~ 3 fi 2 ' 

Further, with uj rn = fl m h 2 and h = Ho/ (100km/ s/Mpc), one writes 



(18) 



(19) 



Similarly, uj ra d is defined for radiation but we consider its contribution to be negligible at late times. 

In terms of the logarithmic variables, the generalized Friedmann equations (equations (|14l) (|15[) (|16[1 ) may be written 



Model DFT-A: 



Model DFT-B: 



(e 2 " - HDifiG 2 + G(u) 2 )) - \VG~o\({G 2 + G(u) 2 ) In ( ° l + ^ {u) 

-1152/t 8 e 8tI (4u , (u' + 1) + «")) = 0, 
(e 2 --Hlm 2 {Gl + G{u) 2 ) 2 )-G{uf 
+576/* 8 e 8u G 2 l ((7u' - !)(«' + 1) + 2u") 



Model DFT-C: (e 



2u tt2 \ 



= 0, 



G A / -576/i 8 e 8ll (4wV + 1)) + u") 



3/t 



■(■ 



Gn cosh 2 



PS?) 



(20) 
(21) 

(22) 
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where H 2 = (Cj r e- iN + u m e~ 3N ) + ^ and G(u) = ^ 4 e 4 "24(u' + 1). 

Unlike equations (TT4")) , (|15p , and (|16p , we found that the equations written in terms of the logarithmic variables allow 
stable numerical integrations over redshift ranges from z = 5 to z — 0. As discussed in (l9| . it is necessary to perform 
the integration forward in time (backward in the redshift) with initial conditions provided by some approximate 
solutions at high redshift. We verified that at higher redshifts (z > 5), the approximate solutions provide an excellent 



fit to the numerical solution of the ODEs with a relative difference 



that is better than 0.1%. We use 



these approximations in order to find the numerical initial conditions for u and u' at z — 5 and then start a numerical 
integration of the ODEs down to the supernovae locations. Proceeding in this way, we obtained very stable programs 
to derive various Hubble plots for the HOG models. The approximate solutions for the models of this section can be 
derived as 



'GqX 



Model DFT-A: H 2 = + + g( , )a) 



((G 2 +G(?i) 2 )ln( 



r,\2, 



Gj + Gj u) 

r 2 



Model DFT-B: H 2 = jjfH 



Model DFT-C: H 2 = tfH% 



1152/ t 8 e 8fi (4u'(w' + l)+u")^ 

G{uf - 576fL s e m G 2 ((7u' - l)(u' + 1) + 2u") 



i{Gl + G{u) 2 Y 
fG^X /-576/2 8 e 8il (4u'(u' + 1)) + u") 



(23) 
(24) 



G(u) 
+ — ^tanh 
Go 
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G 2 cosh 2 (<f> 



-In 



cosh 



Go 



(25) 



where G(u) = /t 4 e 4tl 24(u' + 1) and the source u is used from earlier. 

The generalized Friedmann equations above along with their high redshift approximations are in a form ready for 
numerical integration and comparisons to observational data in §4. 



3.2. f(G) models proposed by Zhou, Copeland and Saffin 



The authors of 22j performed a thorough phase space analysis of f(G) models and analyzed some specific models 
that satisfy some cosmological viability conditions. Following [23[ for f(R) models, the authors of [22[ expressed the 
viability conditions as constraints on the derivatives of the function f(G). We consider here some of their models as: 



Model ZCS-A 
Model ZCS-B 
Model ZCS-C 



f(G)^aVG + l3^G, 
/(G) = a(G 3 / 4 -/?) 2 / 3 , 
/(G) = ^exp(/3/G). 



where a and f3 are constants. The equations of motion in a flat FLRW spacetime follow as: 
Model ZCS-A: 



■AH 2 = -~(\aVGH{H 3 ~ \{H' 2 H + H" H 2 )) 



~(3VGH(^H' 2 H + H 3 + H 2 H' - ^{H /2 H + H"H 2 ))j / (h'H + H 2 



(26) 
(27) 
(28) 



a3N 



a 4N 



(29) 



Model ZCS-B: 



3H 2 = -a[2(H 2 (H'H + H 2 )) 3 ^ (3 3 H(H' + H) - 576H 12 {H' + H)2A 1 ' A 

--H 4 (9H' 2 + 18H'H + WH 2 )(H 2 (H'H + H 2 )) 1/2 2A 3/i 
+12j3H 7 (5H' 3 + 15H + 18H'H 2 + 8H 2 )(H 2 (H'H + i? 2 )) 1/4 24 1/2 
+ (288H 9 (H' + iT) 2 24 1/4 + ^H 3 fi 2 (H 2 {H 1 H + H 2 ))2A 3/i 



■5 



-18H 6 P(H' + H)(H 2 (H'H + i? 2 )) 1 / 4 24 1 / 2 ) (H' 2 H + H"H 2 )^j / 
{{H'H + H 2 ){2A 3 ' 4 {H 2 {H'H + H 2 )) 3 ' 4 - p) 7 ' 3 {H 2 {H'H + if 2 )) 3 / 4 ) 

3/v 1 tin y 6V > 



Model ZCS-C: 



3H 2 = -aexp ( 24g2/(j ^ g + g2) ) ((28817 1 + ff) 2 + 2/3 2 ff 2 ff'(ff' + 2ff) 

+ (H' 2 H + H" H 2 ){2AH 4 /3{H' + H) + /3 2 H - 1UH 7 (H' + H) 2 )){2AH 3 H' + 24H 4 ) 1/2 
+24H^(3(H 3 + 9H' 2 H + 7H'H 2 + 3H' 3 )^ j (h76H\H'H + H 2 ) 4 ^ 

Again, in order to compare these models to cosmological data, we must write the generalized Freidmann equation in 
a numerically stable integrable form using logarithmic variables yielding: 
Model ZCS-A: 



(e 2u - e 2 ")(48AV + l) 2 ) - Aa^G{u){u" + 2{u' + \){u' - 1)) 

-3^(u"-4(u' + l)) = 0, (32) 



Model ZCS-B: 



Model ZCS-C: 



(e 2u - e 2i )6/tV + l)(G(w)(G(u) 3/4 - p) 7/3 ) - (G(u) 13/ V - 1) 
+G{uf /2 p{u' + 4) - G{u) 7/4 p 2 {Au' + 5) + 2G(u)/3 3 {u' + 1)) 

-6fi 4 e 4u {2G(uf' 4 - 3G{uf/ 2 p + G(uf/ 4 p 2 )u" = 0, (33) 

(e 2u - e 2i )(3A 2 K + l) 2 G(w) 2 ) + aty/Gfuj exp(/3/G(u)) (( - ^G(u)V - 1) 

+G(u)0(5u' + 1) + 4/3V) (u' + 1) + (> - ic( ? i) 2 + G(u)/?) u") = 0, (34) 

We also need the approximate solution for these models at high redshift providing initial conditions for numerical 
integrations as discussed in previous section. The approximate solutions at high redshift for these models are: 
Model ZCS-A: 



_ 4a^W)(u" + 2(u' + l)(u' 1)) + 3(3^/W)(u" 4(u> + 1)) 

^ ~ /ie + 48(£' + l) 2 ' {Sb) 

Model ZCS-B: 

if 2 = A 2 e 2 " + ((G(u) 13 / 4 (?i'-l) + G(u) 5 / 2 /3( t i' + 4)-G(u) 7 / 4 /3 2 (4u' + 5) 
+2G(w)/3 3 (w + 1)) + 6^ 4 e 4fi (2G(G) 9/4 - 3G(u) 3/2 /3 

+G(£) 3 / 4 /3 V') /(q{l{G{uY' 4 - P) 7 ' 3 G{u){u' + 1)) , (36) 

Model ZCS-C: 

H 2 = fi 2 e 2Z + ay/G(u)exp(P/G(u)) ( ( - ]^G{u) 2 {u' - 1) + G{u)/3(5u' + 1) + 4/3 2 u') (u' + 1) 

+ (/? 2 - \g(u) 2 + G(u)p)u") /(3AV + l) 2 G(w) 2 ), (37) 
where again G(u) = fi 4 e 4u 2A(u' + 1) and the source u is as defined earlier. 
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FIG. 1: LEFT: 2D joint contour plots for DFT-A with parameters Q m , a, A, p., H for Union, WMAP5 and SDSS LRG(BAO) 
data sets, where the inner and outer loops are 68% and 95%, respectively. RIGHT:2D joint contour plots for DFT-A with 
parameters fi m , a, A, Go, Ho for Union, WMAP5 and SDSS LRG(BAO) data sets, where the inner and outer loops are 68% 
and 95%, respectively. 



3.3. f(G) Models proposed by Uddin, Lidsey and Tavakol 

The models presented by [24| are similar to the models presented by [22j of the previous section, i.e. /(G) = 
ay/G + p-\/G, with = and a = y/£a, although 24] performed a different and independent analysis. In agreement 
with the analogous scalar field models, the power-law form of the /(G) model as /(G) = ±2v / aG was found to 
have stable scaling solutions. The authors studied the equation of state parameter for these models as it evolved 
through radiation, matter and accelerating epochs. The behavior of the energy densities were also discussed there. 
By investigating stability for both vacuum and non-vacuum solutions, it was recognized that those models suffered 
from a singular point at transition. Our results for these models can be found with results for the ZCS models in §5 
below. 



4. COSMOLOGICAL DISTANCES CONSTRAINTS FROM SNE IA, CMB SURFACE, AND BAO 

In this section, we describe the three cosmological observations used to constrain the models described above. One 
of the first compelling evidences for cosmic acceleration came from the Supernovae type la (SNe la) observations and 
we use here the distance modulus as a function of the redshift z given by 

fi(z) = m- M = 5\ogD L +25 (38) 

containing the magnitude-redshift function, rh(z) and a nuisance parameter, M which is degenerate with the Hubble 
parameter, Hq. The luminosity distance, Dl in units of Mpc for HOG models is given by 

D L (z) = (l + z)£ j^Kdz', (39) 

with H(z') as the solution to the non-linear differential equations given above for the higher order gravity models 
(recall that for the logarithmic variables defined in the previous section the redshift reads as z — e~ N — 1.) We use the 
Union set of supernovae which was compiled as an attempt to gather the best of the best supernovae from different 
surveys, including Supernovae Legacy Survey, ESSENCE Survey, HST, and other older sets [25|]. After selection cuts, 
the 414 SNe la are reduced to 307. The fitting of the SNe la uses a standard minimization method, x 2 , 

i=l 
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FIG. 2: LEFT: 2D joint contour plots for DFT-B with parameters Q m , a, X, G , H for Union, WMAP5 and SDSS LRG(BAO) 
data sets, where the inner and outer loops are 68% and 95%, respectively. RIGHT:2D joint contour plots for DFT-C with 
parameters Q m , a, X, Go, Ho for Union, WMAP5 and SDSS LRG(BAO) data sets, where the inner and outer loops are 68% 
and 95%, respectively. 



where a is the magnitude uncertainty and i the number of data points compared. We break degeneracies in the 
parameter space by using the value of the Hubble parameter as measured by the Hubble Key Project (HST)|26], 
H = 72 ± 8 km/s/Mpc. 

Next, we consider the distance to the CMB last scattering surface and, following [27|], we define three fitting 
parameters for comparison to the WMAP5 data using the shift parameter, R, see for example (28j . 

R(z*) = \/n^H a {l + z*)D A (z*), (41) 
with the redshift, z* for the surface of last scattering, see for example [29j |. 

z* = 1048[l + 0.00124(^f,/i 2 )-°- 738 ][l+5i(rj m /i 2 ) S2 ], (42) 

where 

0.0783(^ 2 )- - 238 
91 l + 39.5(fW) - 763 ' 



and 



0.560 

92 i + 2i.i(n 6 ft a ) 1 - 81 ' 

and third, the acoustic scale, l a , is, see for example [3(1 l3lj] . 

l a = (l + Z *f-^, (45) 

with the proper angular diameter distance, Da{z) = Dl{z)/(1 + z) 2 and the comoving sound horizon, r g (z*), see for 
example [321 ] . 

r ( z ) — / (46) 

VZJo a 2 H{a)^/\T(^hl^Q^' 

with f2 7 = 2.469 x 10 ~ 5 /i~ 2 for T CTO & = 2.725-ftT. Together the parameters = (R,l a ,z^) are used to fit XwmAP = 
AxiCov~ 1 (xiXj)Axj with Ax, — Xi — x° bs and Cov~ 1 (xiXj) is the inverse covariance matrix for the parameters. 

Thirdly, we use BAO constraints following [13, 132rl34l | and defining the ratio of the sound horizon, r s (zd) to the 
effective distance, Dy as a fit for SDSS by 

2 _( r s (z d )/D v {z = 0-2) -0.198 ^ ( r s (z d )/D v (z = 0.35) - 0.1094 ^ 
Xbao ~ V O0058 J + V O0033 J ' (4?) 



8 




0.22^^^^^^^H 0.221 

-0.02 -0.01 0.01 -0.06 -0.04 -0.02 





0.22^^^^^^^H 

1.2 1.3 1.4 1.5 70 72 74 

u H. 




-2-1 1 



-0.4 -0.2 





1 .25 1 .3 1 .35 



70 71 72 73 




FIG. 3: LEFT: 2D joint contour plots for ZCS-A with parameters Q. m , a, /3, H for Union, WMAP5 and SDSS LRG(BAO) 
data sets, where the inner and outer loops are 68% and 95%, respectively. CENTER:2D joint contour plots for ZCS-B with 
parameters fi m , a, f3, Ho for Union, WMAP5 and SDSS LRG(BAO) data sets, where the inner and outer loops are 68% and 
95%, respectively. RIGHT:2D joint contour plots for ZCS-C with parameters fi m , a, 0, H for Union, WMAP5 and SDSS 
LRG(BAO) data sets, where the inner and outer loops are 68% and 95%, respectively. 



with, see for example [33I ]. 



and the rcdshift, Zd as, see for example 



D v {z) = {D\{z){l + zf 1 ^) 1/ \ 



where 



and 



i29i(q»fo 2 ) - 251 2b 



bi = 0.313(r> m / l 2 )-°- 419 [l + 0.607(O m /i 2 ) Ub74 ] 



b 2 = 0.238{Q m h 2 ) 



2\0.223 



(48) 



(49) 



(50) 



(51) 



We use a modified and extended version of the publicly available package, CosmoMC [36[ to perform a Monte 
Carlo Markov Chain analysis (MCMC's). The MCMC's are used to compute the likelihoods for the parameters in 
the model. This method randomly chooses values for the above parameters and, based on the x 2 obtained, either 
accepts or rejects the set of parameters via the Metropolis-Hastings algorithm. When a set of parameters is accepted 
it is added to the chain and forms a new starting point for the next step. The process is repeated until the specified 
convergence is reached. We found that an elaborate and essential step in our analysis was to numerically integrate the 
stiff ODEs in order to derive the Hubble expansion rate for a wide range of redshift and integrate it to the CosmoMC 
package. For this task, we made a necessary change of variable in our integration and also used a good approximation 
at higher redshifts in order to obtain initial values for our codes to perform stable integrations from z = 5 down to 
z = 0. We used the framework developed to compare the models described in the previous sections to the three sets 
of observations and the results are given in the next section. 



RESULTS AND CONCLUSION 



Our results are summarized in figures 1-4 and in table I. In figure 1, we present the results for model DFT-A. We 
first (left part of figure 1) varied fi m , a, A, fi, and Hq leaving Go ~ /i 4 as a derived parameter, and we obtained a 
normalized x 2 /do/ = 1.0754 (x 2 — 329.1). On the right part of figure 1 for model DFT-A, we varied Q, m , a, A, Go, 
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FIG. 4: Comparisons of different SNe la plots with Union data and different values of f(G) models. LEFT: Comparison of 
DFT-A, DFT-B, DFT-C Models and LCDM. RIGHT:Comparison of ZCS-A, ZCS-B, ZCS-C Models and LCDM. 
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TABLE I: x 2 ) X 2 /dof, and best-fit parameters for f(G) higher order gravity models using observational constraints from 
supernovae, HST, CMB surface and BAO. For models DFT-A, DFT-B, DFT-C the parameters are Q m , H , (I, Go, a, and A. 
For models ZCS-A, ZCS-B, ZCS-C the parameters are fl m , Ho, fl, a, and /3. 



and H leaving fx as a derived parameter, with a normalized x 2 /dof = 1.0751 (x 2 = 329.0). Our best fit values are 
given in table 1 and we find an fi m around 0.25 in both cases. We find best values for the model parameters a and A 
that have a product aA of the order unity as required, see [20T | . We find the best fit for fx that is of the order of twice 
the normalized Hubble parameter h, as it should 0, [2l|, and the best fit parameter Go of the order of h 4 [2(|. Next, 
we show in figure 2 the results for models DFT-B and DFT-C with similar results and normalized x 2 /dof — 1.0764 
(x 2 = 329.4) and x 2 /dof = 1.0748 (x 2 = 328.9) respectively. Finally, we show in figure 3, results for the ZCS models 
varying the parameters a, /3, (x, f2 m and H . 

The best fit parameters are given in table 1. The first and second models have best fit values similar to the DFT 
models while the third model has slightly higher values for Q, m and Hq. The best fit x 2 ' s f° r f{G) models are close 
to the x'sn+bao+cmb ~ 314.6 (x 2 /dof — 1.018) that we obtained for the LCDM concordance model. In view of the 
of the possible systematic uncertainties in the supernova data, it is not clear that the difference between the x 2 s is 
significant. It is worth noting that the parameter space that we found for the DFT models is also all contained within 
the parameter space found by [37j for the models to be compatible with solar system constraints. In other words, the 
DFT models analyzed in this paper pass physical acceptability conditions [2(J [H| , solar system tests [37f and here pass 
constraints from supernova, BAO rulers and distance to CMB last scattering surface. Recently, Ref. [39| pointed out 
in a preliminary work that matter perturbations in Gauss-Bonnet models exhibit some instabilities during the matter 
era, and that for the growth to be compatible with observations the deviations from general relativity have to be very 
small. This point needs further investigations using perturbation studies in /(G) models. In view of the success of 
/(G) models with solar system tests and cosmological distances constraints, we conclude that these models need to 
be subjected, in future projects, to full large scale structure constraints such as galaxy clustering and gravitational 
lcnsing, as well as the full CMB analysis. 
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